clearvars;

% =========================================================================
% 超声和光声的SIMULATION，光声仿真时的初始压力分布跟超声仿真有何不同
% =========================================================================

% create the computational grid
PML_size = 20;          % size of the PML in grid points
Nx = 440 - 2*PML_size;  % number of grid points in the x (row) direction
Ny = 240 - 2*PML_size;  % number of grid points in the y (column) direction
dx = 0.1e-3;            % grid point spacing in the x direction [m]
dy = 0.1e-3;            % grid point spacing in the y direction [m]
kgrid = kWaveGrid(Nx, dx, Ny, dy);

%% 10MHZ 探头参数==================%
%% PA&US shared reconstruction parameters =========%
params.Nelements = 64;
params.fs = 40e6;   % [Hz]
params.fc = 10e6;  % [Hz]
params.bandwidth = 90;    % [percent]
params.pitch = 0.15e-3;  % pitch == spacing [mm]
params.width = 0.127e-3;  %  [mm]
params.radius = inf;   % 线阵
params.c = 1500;  %[m/s]
params.fnumber = [];
xlen = params.pitch * 63 + params.width;
x = (-xlen/2):0.025e-3:(xlen/2);
z = 0:0.025e-3:22.5e-3;  % exp 60mm = 60e-3
[X, Z] = meshgrid(x, z); % 生成网格点坐标[mm]

% define the properties of the propagation medium
% Define the acoustic properties
%globe
medium.sound_speed = 1500*ones(Nx, Ny);    % [m/s]
medium.density = 1000*ones(Nx, Ny);        % [kg/m^3]
%bone
medium.sound_speed(find(bone==1)) = 2500;   % [m/s]
medium.density(find(bone==1)) = 1300;   % [kg/m^3]
%Scatter
medium.sound_speed(find(discs==1)) = 2500;   % [m/s]
medium.density(find(discs==1)) = 1300;   % [kg/m^3]

% create initial pressure distribution using makeDisc

%%随机点状结构
% disc_magnitude = 5; % [Pa]
% disc_x_pos = 60;    % [grid points] 可以设置为随机位置
% disc_y_pos = 60;  	% [grid points]
% disc_radius = 5;    % [grid points]
% disc_2 = disc_magnitude * makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);
% 
% disc_x_pos = 30;    % [grid points]
% disc_y_pos = 30; 	% [grid points]
% disc_radius = 8;    % [grid points]
% disc_1 = disc_magnitude * makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);
% 
% source.p0 = disc_1 + disc_2;

% % 板状结构
% disc_magnitude = 20;
% angle = pi;
% length = 40;
% thickness = 10;
% start_point = [49, 30];
% source.p0 = zeros(Nx, Ny);
% for i = 1 : thickness
%     start_point = [49+i, 30];
%     line = disc_magnitude * makeLine(Nx, Ny, start_point, angle, length);
%     source.p0 = source.p0 + line;
% end

% % 弧状结构
% disc_magnitude = 20;
% disc_x_pos = 60;    % [grid points] 可以设置为随机位置
% disc_y_pos = 60;  	% [grid points]
% disc_radius1 = 10;    % [grid points]
% disc_1 = disc_magnitude * makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius1);
% magnitude2 = 10;
% disc_x_pos2 = 60;    % [grid points]
% disc_y_pos2 = 70; 	% [grid points]
% disc_radius2 = 15;    % [grid points]
% disc_2 = magnitude2 * makeDisc(Nx, Ny, disc_x_pos2, disc_y_pos2, disc_radius2);
% disc = disc_1 + disc_2;
% disc(disc ~= disc_magnitude) = 0;  % 使用逻辑索引将不等于 disc_magnitude 的元素置为 0
% source.p0 = disc;

% % 血管结构  STARE dataset
% load the initial pressure distribution from an image and scale the
% magnitude
disc_magnitude = 20;
p0 = disc_magnitude * loadImage('./POINTS_Dataset/points_initial_pressure_0001.png');
% resize the image to match the size of the computational grid and assign to the source input structure
source.p0 = resize(p0, [Nx, Ny]);


% smooth the initial pressure distribution and restore the magnitude
source.p0 = smooth(source.p0, true);

% define a binary line sensor
sensor.mask = zeros(Nx, Ny);
sensor.mask(1, :) = 1;   %  设置成64
% 设置探测器的响应的频率带宽
%sensor.frequency_response

% create the time array
kgrid.makeTime(medium.sound_speed);

% set the input arguements: force the PML to be outside the computational
% grid; switch off p0 smoothing within kspaceFirstOrder2D
input_args = {'PMLInside', false, 'PMLSize', PML_size, 'PlotPML', false, 'Smooth', false};

% run the simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});

% 增加有限带宽滤波，也可以通过设置sensor.frequency_response实现
% filter the sensor data using a Gaussian filter
Fs = 1/kgrid.dt;        % [Hz]
center_freq = 2.5e6;      % [Hz]
bandwidth = 100;        % [%]
sensor_data_gaussian = gaussianFilter(sensor_data, Fs, center_freq, bandwidth);

% reconstruct the initial pressure
p_xy = kspaceLineRecon(sensor_data.', dy, kgrid.dt, medium.sound_speed, 'Plot', true, 'PosCond', true, 'Interp', '*linear');
%p_xy = kspaceLineRecon(sensor_data_gaussian.', dy, kgrid.dt, medium.sound_speed, 'Plot', true, 'PosCond', true, 'Interp', '*linear');

% define a second k-space grid using the dimensions of p_xy
[Nx_recon, Ny_recon] = size(p_xy);
kgrid_recon = kWaveGrid(Nx_recon, kgrid.dt * medium.sound_speed, Ny_recon, dy);

% resample p_xy to be the same size as source.p0
p_xy_rs = interp2(kgrid_recon.y, kgrid_recon.x - min(kgrid_recon.x(:)), p_xy, kgrid.y, kgrid.x - min(kgrid.x(:)));

% =========================================================================
% SAVE RESULTS，需要保存ground-truth图像，采集到的超声信号，光声信号，粗重建得到的超声图像，光声图像
% =========================================================================

% 保存 ground-truth图像
% 定义保存路径和文件名
filename = 'ground_truth_initial_pressure.png';
% 将 source.p0 的数据范围缩放到 [0, 255] 以便于保存为图像
scaled_p0 = mat2gray(source.p0);
% 使用 imwrite 函数保存图像
imwrite(uint8(scaled_p0 * 255), filename);

% 保存 采集到的光声信号
save('simulation_photoacoustic_results.mat', 'sensor_data');
%save('simulation_photoacoustic_results.mat', 'sensor_data_gaussian');

% 保存 粗重建得到的光声图像
filename = 'reconstructed_photoacoustic.png';
% 将 source.p0 的数据范围缩放到 [0, 255] 以便于保存为图像
scaled_rs = mat2gray(p_xy_rs);
% 使用 imwrite 函数保存图像
imwrite(uint8(scaled_rs * 255), filename);

% =========================================================================
% VISUALISATION
% =========================================================================

% plot the initial pressure and sensor distribution
figure;
imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, source.p0 + sensor.mask * disc_magnitude, [-disc_magnitude, disc_magnitude]);
colormap(getColorMap);
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;
colorbar;
title('Initial Pressure and Sensor Distribution');
scaleFig(1, 0.65);

% plot the simulated sensor data
figure;
imagesc(sensor_data, [-1, 1]);
colormap(getColorMap);
ylabel('Sensor Position');
xlabel('Time Step');
colorbar;
title('Simulated Sensor Data');

% plot the reconstructed initial pressure 
figure;
%imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, p_xy_rs, [-disc_magnitude, disc_magnitude]);
imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, p_xy_rs, [-disc_magnitude, disc_magnitude]);
colormap(getColorMap);
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;
colorbar;
title('Reconstructed Initial Pressure'); 
scaleFig(1, 0.65);

% plot a profile for comparison
figure;
plot(kgrid.y_vec * 1e3, source.p0(disc_x_pos, :), 'k-', kgrid.y_vec * 1e3, p_xy_rs(disc_x_pos, :), 'r--');
xlabel('y-position [mm]');
ylabel('Pressure');
legend('Initial Pressure', 'Reconstructed Pressure');
axis tight;
set(gca, 'YLim', [0, 5.1]);
title('Pressure Profile Comparison'); 